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SUMMARY 

We study rotating thermal convection in spherical shells as prototype for flow in the cores of 
terrestrial planets, gas planets or in stars. We base our analysis on a set of about 450 direct nu¬ 
merical simulations of the (magneto)hydrodynamic equations under the Boussinesq approxi¬ 
mation. The Ekman number ranges from 10“^ to 10“^. The supercriticality of the convection 
reaches about 1000 in some models. Four sets of simulations are considered: non-magnetic 
simulations and dynamo simulations with either free-slip or no-slip flow boundary conditions. 
The non-magnetic setup with free-slip boundaries generates the strongest zonal flows. Both 
non-magnetic simulations with no-slip flow boundary conditions and self-consistent dynamos 
with free-slip boundaries have drastically reduced zonal-flows. Suppression of shear leads 
to a substantial gain in heat-transfer efficiency, increasing by a factor of 3 in some cases. 
Such efficiency enhancement occurs as long as the convection is significantly influenced by 
rotation. At higher convective driving the heat-transfer efficiency tends towards that of the 
classical non-rotating Rayleigh-Benard system. Analysis of the latitudinal distribution of heat 
flow at the outer boundary reveals that the shear is most effective at suppressing heat-transfer 
in the equatorial regions. Simulations with convection zones of different thickness show that 
the zonal flows becomes less energetic in thicker shells, and, therefore, their effect on heat- 
transfer efficiency decreases. Furthermore, we explore the influence of the magnetic field on 
the non-zonal flow components of the convection. For this we compare the heat-transfer ef¬ 
ficiency of no-slip non-magnetic cases with that of the no-slip dynamo simulations. We find 
that dii E = 10“^ magnetic field significantly affects the convection and a maximum gain of 
about 30% (as compared to the non-magnetic case) in heat-transfer efficiency is obtained for 
an Elsasser number of about 3. Our analysis motivates us to speculate that convection in the 
polar regions in dynamos at E" = 10“^ is probably in a ‘magnetostrophic’ regime. 
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1 INTRODUCTION 

The magnetic field we see on the Earth, on other planets and stars, 
is thought to be generated via a dynamo mechanism which is pow¬ 
ered by the convective motions in the interior of these objects (e.g. 
Brandenburg and Subramanian 2005; Donati and Landstreet 2009; 
Wicht and Tilgner 2010; Jones 2011; Roberts and King 2013). 
In terrestrial planets like the Earth, the convection is driven by a 
combination of compositional and thermal buoyancy, while for gas 
planets and stars convection of thermal origin is a more efficient 
means of heat transport than conduction or radiation. 
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Due to the advances in the computing power numerical sim¬ 
ulation of the equations that govern convective flow has become 
an important tool for enhancing our understanding of convection in 
natural objects. A lot of effort has been dedicated to studying ther¬ 
mal convection in the simplified setup of a plane layer where the 
bottom boundary is hotter than the top boundary (e.g. Ahlers et al 
2009). This idealized form of convection is referred to as Rayleigh- 
Benard convection (RBC). Many studies have focused on the heat- 
transfer efficiency of the convection as a function of the various 
physical parameters of the system. An understanding of such a re¬ 
lationship allows us to model the thermal evolution of planets and 
stars (e.g. Breuer et al. 2010; Kippenhahn et al. 2012), and can help 
to elucidate the nature of the dynamo operating in such systems. 

Earth, other planets, and stars rotate and the associated Corio¬ 
lis forces affect the properties of the convection occurring in their 
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interior. The Proudman-Taylor theorem (Proudman 1916; Taylor 
1923) demonstrates that in rotating systems, Coriolis forces tend to 
suppress motion in the direction along the rotation axis, rendering 
the flow nearly two-dimensional. If the RBC setup described above 
rotates along the vertical direction, convection occurs in columnar 
cells aligned with the rotation axis. The vertical component of the 
flow, necessary for heat transfer, will be signiflcantly reduced. 

The presence of magnetic held is ubiquitous in planets and 
stars. It is well known that imposing a magnetic held suppresses 
the convective motions via the associated Lorentz forces in a typi¬ 
cal RBC setup (Chandrasekhar 1961). Magnetic quenching can re¬ 
duce the heat-transport efficiency of the convecting flow. The effect 
of the dynamo-generated magnetic held on the convective flow is 
more complex than that of a uniform imposed held and depends 
on various factors, for example, the rotation rate of the system, the 
degree of density stratiflcation in the convecting layer, etc. 

The interplay of the effect of rotation and magnetic held on 
the convection is of particular importance as it relates to the situ¬ 
ation in planets like the Earth which rotate rapidly and possess a 
strong dynamo-generated magnetic held (Roberts and King 2013). 
Acting alone, both rotation and magnetic held quench the con¬ 
vection. However, under some conditions, the presence of mag¬ 
netic held in a rotating system can enhance the heat-transport ef- 
flciency as compared to the corresponding rotating setup without 
magnetic held (Chandrasekhar 1961; Stevenson 1979). This ‘mag- 
netostrophic’ regime, in which Lorentz force and Coriolis force are 
(aside from pressure forces) dominant in the force balance of the 
flow, is thought to be of fundamental importance for the Earth and 
for other rapidly rotating planets (Jones 2011; Roberts and King 
2013; Davidson 2013). 

Geometry can also play a crucial role in determining the be¬ 
haviour of convecting fluids. It is well known that convection in 
rotating spherical shells can self-consistently excite strong axisym- 
metric azimuthal (zonal) flows (Gilman 1977), with free-slip me¬ 
chanical boundary conditions being especially suitable for per¬ 
mitting such flows (Gilman 1978). The excitation of zonal flows 
in rotating spherical shells by Reynolds stresses is thought to 
be due to the effect of the curved spherical boundaries on the 
convection (Busse and Hood 1982). In non-rotating two dimen¬ 
sional RBC with horizontally periodic boundaries similar shear 
flow can be excited via the so-called ‘tilting’ instability of convec¬ 
tion rolls (Thompson 1970; Eitzgerald and Earrell 2014). Although, 
the non-rotating three-dimensional (3D) RBC with periodic bound¬ 
aries does not generate strong zonal flows. In such systems, even if 
localized tilting of convection cells excite some zonal flow, simi¬ 
lar processes in neighbouring locations may suppresses it since the 
convective rolls might be randomly oriented. However, such 3D se¬ 
tups can promote strong zonal flows if they rotate along a horizontal 
direction (Hardenberg et al 2015). Apart from not contributing to 
heat-transfer, such zonal flow is capable of deflecting and shearing 
the heat-transferring flow itself (Terry 2000; Goluskin et al 2014). 

The effects described above highlight that (mag- 
neto)hydrodynamic convection in rotating spherical shells is 
richer in physical phenomena than the classical RBC setup. 
Rotating RBC has been investigated using both theoretical and 
laboratory experiments (see e.g. King et al 2012; Ecke and 
Niemela 2014; Cheng et al 2015, and references therein). One 
of the main results of these studies is the existence of a transition 
regime (see Aurnou et al 2015). Depending on the rotation rate 
and other control parameters of the system, the Coriolis forces can 
either strongly affect the flow and hence the heat-transfer efficiency 
of convection, or other forces (buoyancy or inertia) can dominate 


and the system approaches the non-rotating RBC state. The 
exact relationship describing the dependence of the heat-transfer 
efficiency on the physical parameters of the system is still a matter 
of debate. Rotating RBC with either an imposed magnetic held or 
a self-generated one has also been studied. Stellmach and Hansen 
(2004) reported that applied/self-generated magnetic held could 
substantially change the length scale and vigor of the convection. 
Recently, King and Aurnou (2015) experimentally investigated 
rotating RBC with an imposed magnetic held and also reported 
enhancement of convective efficiency depending on the applied 
magnetic held. 

Rotating and magnetic convection in spherical shells has been 
investigated in the context of geodynamo models. Here, the fo¬ 
cus is usually on characterising the dynamo-generated magnetic 
held. Christensen et al (1999) investigated a set of spherical dy¬ 
namo simulations and explored the effect of various control param¬ 
eters on the strength and morphology of the magnetic held. There 
have been more extensive systematic parameters studies since then 
which have tried to focus on various other relevant aspects. Eor ex¬ 
ample, the effect of the fluid Prandtl number (Simitev and Busse 
2005) and the effect of different velocity and temperature bound¬ 
ary conditions (Christensen 2002; Busse and Simitev 2006; Willis 
et al 2007; Schrinner et al 2012; Davies et al 2013; Yadav et al 
2013b; Dharmaraj et al 2014), on convection and magnetic held 
have been investigated using parameter studies. 

In this study we take a somewhat different approach and in¬ 
vestigate exclusively how convective flow is affected by different 
velocity boundary conditions, control parameters, and the presence 
of dynamo-generated magnetic held. We will build our study on an 
extensive set of numerical simulations which reach parameter val¬ 
ues that were so far hard to access in numerical simulations and 
cover a wide enough parameter space to highlight a few interesting 
and novel trends. The basic setup of the majority of the simulations 
in the data-set is similar to the typical geodynamo models studied 
in the past. 

At this juncture it is appropriate to point out what our study 
does not touch upon. To determine the accurate relationship be¬ 
tween the amount of heat transferred via convective motions and 
the system control parameters is of fundamental importance, but is 
not the aim of this study. In the context of magnetic convection in 
rotating spherical shells, this has been investigated by King et al 
(2010). Here we concentrate on how heat-transfer is affected by 
the factors discussed above: boundary conditions, zonal flow, and 
presence/absence of magnetic held. Eurthermore, recent develop¬ 
ments in the rotating convection theory and simulations (see e.g. 
Julien et al 2012; Cheng et al 2015) propose control parameters 
for asymptotic regimes which are beyond what we can achieve with 
our computational resources. Nonetheless, we expect that at least 
the qualitative trends we And in our simulations are relevant in the 
asymptotic regime in planets and stars. 

We describe the governing equations and the model ingredi¬ 
ents in Sec. 2. In Sec. 3, we analyse purely hydrodynamic con¬ 
vection in rotating spherical shells with different boundary condi¬ 
tions. We then discuss self-consistent dynamo simulations in Sec. 4 
and use two different setups with free-slip boundaries and no-slip 
boundaries to disentangle the effect of magnetic held on zonal and 
non-zonal flow components. We conclude and discuss possible ex¬ 
tensions of our study in Sec. 5. 
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2 MODEL SETUP 

2.1 MHD Equations 


We solve the fundamental (magneto)hydrodynamic (MHD) equa¬ 
tions in a rotating spherical shell containing an incompressible fluid 
of constant density p. The spherical shell’s inner and outer bound¬ 
aries are located at radius Vi and Tq respectively. The ratio Vi / Tq 
deflnes the aspect ratio p of the spherical shell. The spherical shell 
rotates with a frequency Q along the vertical axis z. We work with 
a non-dimensional set of equations where the fundamental units 
are: the shell thickness D — Vo — ri for length scale, the viscous 
diffusion time jv {y is the viscosity) for time scale , yjppXQ. 
(p is the magnetic permeability and A is the magnetic diffusivity) 
for magnetic held scale, and the imposed temperature contrast AT 
between the inner and the outer boundary for temperature scale. 

With these non-dimensional units, the Boussinesq equations 
governing the velocity u, the temperature T, and the magnetic held 
B, are given by 

e(^^ + u-VuJ + 2zxu = -VP+^^g{r)Tr 
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where P is the pressure, and g{r) is the gravity profile which 
is either r/vo or (ro/r)^. In purely hydrodynamical simulations, 
Eqns. 4 and 5 are not relevant and B is set to zero in Eqn. 1. 
The non-dimensional control parameters appearing in the equations 
above are defined as follows: 
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where k is the thermal diffusivity, a is the thermal expansivity, and 
Qo is the gravity at the outer boundary. 


2.2 Boundary Conditions 

We assume the same boundary conditions for both boundaries in 
the majority of the cases. This is done for the sake of simplicity 
and also since we are not modelling a specific planet where such 
a choice might not be appropriate (for instance, in the Earth, in¬ 
ner solid core and outer mantle impose different magnetic bound¬ 
ary conditions). The boundary conditions are also similar to many 
earlier studies (mainly in the geodynamo context) which facilitates 
direct comparison. 

Eor all of the simulations in this study the temperature is fixed 
on both boundaries, with the inner boundary being hotter than the 
outer one by AT. Eor most simulations both boundaries are either 


free-slip or no-slip. Eor simulations with the deepest convection 
zone (?7 = 0.2), the inner flow boundary condition is no-slip while 
the outer one is free-slip (data set was taken from an earlier study). 
The inner core and the space beyond the outer boundary are con¬ 
sidered insulating, i.e. the magnetic held matches a potential held 
on both boundaries for all of the dynamo simulations. However, for 
cases with p = 0.2 the inner core is conducting. 


2.3 Numerical Method 

We use the pseudo-spectral open source^ code MagIC (Wicht 2002) 
to solve the MHD equations (1-5) with the boundary conditions 
mentioned above. MagIC has been successfully tested against the 
community benchmarks, for both Boussinesq (Christensen et al 
2001) and anelastic (Jones et al 2011) convective dynamo sim¬ 
ulations. To ensure that the velocity and the magnetic held are 
solenoidal they are written as a combination of a poloidal and a 
toroidal part as 


u = VxVxCf + Vx Dr, (10) 

B = VxVxPf + Vx Qf. (11) 

The temperature, pressure, and the scalar potentials C,D,P, Q, are 
expanded in terms of the spherical harmonics for the azimuthal 
and the latitudinal direction and as a sum of Chebyshev polyno¬ 
mials for the radial direction. A combined implicit and explicit 
time-advancing strategy is used for time integrating the equations: 
an explicit second-order Adams-Bashforth scheme is used to time- 
advance the non-linear and the Coriolis term and an implicit Crank- 
Nicolson scheme is used for the other terms (Glatzmaier 1984; 
Christensen and Wicht 2007). 


2.4 A parameter study 

In our simulations, the Ekman number ranges from 10~^ to 10“^. 
The fluid Prandtl number is kept constant at unity. In most cases, 
the magnetic Prandtl number was also unity. Three aspect ratios, 
0.2, 0.35, 0.6, are considered to exemplify very deep, moderately 
deep, and moderately thin convection zones. Note that the grav¬ 
ity is changed accordingly to better approximate the distribution of 
matter: g{r) is r/ro in shells with p of 0.2 and 0.35, which may 
represent the case of Earth’s core at present and at some point in 
the past, while it is (ro/r)^ in shells with p of 0.6, which could 
be representative of the convection zone of a solar-type star. The 
data for ?7 = 0.35 is available in tabulated form in the online ver¬ 
sion. More details on the cases with 77 = 0.6 will be published in a 
forthcoming study. 

We have tried to avoid using imposed symmetries in the az¬ 
imuthal direction as much as possible, especially in the dynamo 
simulations. However, to expedite the initial transient stages of 
many demanding simulations we used either 2-fold or 4-fold sym¬ 
metries in the beginning. The results from these simulations were 
then tested in full sphere calculations. The resolution requirements 
were checked by analysing the resulting kinetic and magnetic en¬ 
ergy spectra (at least 3 decades of decay), by ensuring that the 
time-averaged total heat flux on the inner and the outer boundary 
match, and by matching the input buoyancy power and the dissi¬ 
pated power. 
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Figure 1. Nusselt number Nu as a function of the Rayleigh number Ra 
for hydrodynamic simulations with free-slip boundary conditions, shown 
by filled data points, and hydrodynamic simulations with no-slip boundary 
conditions, shown by empty data points. The color/shape represents dif¬ 
ferent Ekman numbers. The two gray-shaded boxes encapsulate the regime 
where zonal-fiow bistability exists for free-slip cases. The data are for spher¬ 
ical shells with 77 = 0.35. 


5 in Cheng et al. (2015). The cause for the low efficiency of heat 
transfer over some range of Rayleigh numbers is that near the on¬ 
set of convection, columnar cells are arrayed in the vicinity of the 
tangent cylinder (an imaginary cylinder aligned with the rotation 
axis, around the inner non-convective region) and do not fill the 
entire shell volume (Jones et al 2000; Dormy et al 2004). Hence 
conduction must transport the heat over a large part of the distance 
between inner and outer boundary. 

Regime Ri is gradually taken over by a different scaling be¬ 
haviour of the Nusselt number where it rises strongly as Ra in¬ 
creases. This regime is marked with 'R 2 ’. The extent of R 2 is the 
largest in E=10~^ data. A best-fit power law of form Nu oc Ra^ 
gives a ^ 1.23 for data with Ra ^ 2 x 10® and E=10~^. This 
is close to the scaling law Nu oc Ra^^^ proposed by King et al 
(2010) for heat-transfer in rotating spherical-shell convection. As 
we increase Ra further another transition to a different regime oc¬ 
curs where the scaling behaviour of the Nusselt number becomes 
less steep, indicated with marker 'Rs. In this regime, the iner¬ 
tial forces start to dominate over the Coriolis forces. For data with 
Ra ^ 2 X 10® and E=10~^, the power law exponent a is about 
0.35, steeper than the value 2/7 (e.g. King etal 2010) and 1/3 (e.g. 
Ahlers et al 2009). This suggests that our R 3 regime for E=10~^ 
is not asymptotic but rather a ‘bridging’ one; for higher Ra, the 
slope should change to values expected in non-rotating RBC, i.e. 
1/3, or even to 1/2 for very high Ra (Ahlers et al 2009; Gastine 
et al 2015). 


3 HYDRODYNAMIC SIMULATIONS 


3.1.1 Flow structure 


3.1 Models with free-slip boundaries 


The convective heat-transfer is usually measured in terms of the 
ratio of the total heat transported and the heat conducted when con¬ 
vection is absent. This quantity is called the Nusselt number whose 
general definition is 


Nu{r, t) 


In, n~’\ 1 

{UrT)0^(f) 

_ t dTc 
Pr dr 


( 12 ) 


where the angular brackets denote averaging and the subscipts rep¬ 
resent the averaging directions, and dTjdr and dTddr are the to¬ 
tal and conductive temperature gradients, respectively. In the fol¬ 
lowing discussion we will use the time-averaged value of on 
the outer boundary of the spherical shell. Noting that the conduc¬ 
tive temperature gradient on the outer boundary of the shell is — 77 , 
the outer boundary Nusselt number becomes 


Nu = --^ENt(sir = ro). (13) 

77 or 

We begin by considering the heat-transfer behaviour of hydro- 
dynamic (HD) convection with free-slip boundaries in a spherical 
shell with aspect ratio 0.35. The classical Nu vs. Ra behavior is 
shown in Fig. 1 for these cases with filled data points. For A=10“® 
and 10“^ roughly three regimes can be noticed. Just after the onset 
of convection (i.e. where Nu starts becoming larger than 1 ), Nu 
rises very slowly as the supercriticality of the convection increases. 
Marker 'Ri highlights this regime at different Ekman numbers. 
Here the effect of the inertial forces are rather small. This is usu¬ 
ally referred to as the weakly non-linear regime where a scaling of 
type Nu — 1 (X {Ra — Rac)/Rac, where Rac is the Rayleigh 
number at which the convection starts, holds (e.g. Gillet and Jones 
2006; Jones 2007). This regime is peculiar to spherical shells and 
does not exist in plane layer rotating RBC; see, for example. Figure 


In Fig. 2(a, b) we display the axial vorticity (V x u')^ in the equa¬ 
torial plane (as seen from a northern viewpoint) and azimuthally 
averaged zonal flow in a meridional plane for two exemplary 
non-magnetic cases with free-slip mechanical boundaries. Note 
that we have used the non-zonal velocity u' = u — this 

is done to isolate the vorticity contribution of the convective flow 
component. In rotating spherical shells, where convection is signif¬ 
icantly influenced by the Coriolis forces, the convection takes the 
form of alternating cyclonic and anti-cyclonic vortical columns (for 
a review, see Jones 2007; Aurnou et al 2015). Equatorial mani¬ 
festations of such vortical columns are clearly visible in the axial 
vorticity plot in Eig. 2a which shows adjacent regions of positive 
and negative vorticity. In this regime, the columns roughly main¬ 
tain their identity in ^-direction until they hit the outer boundary. 
The prograde tilt (i.e. clockwise spiralling morphology) of these 
columns is also apparent in the vorticity plot. 

In Eig. 2a, the zonal flow profile is almost invariant along 
the rotation axis and the flow near the equator is rotating eastward 
(prograde). Such kind of zonal winds have been investigated exten¬ 
sively in the past, both analytically (e.g. Zhang 1992; Busse 1994) 
and numerically (Gilman 1977; Miesch et al 2000; Grote et al. 
2000; Aurnou and Olson 2001; Christensen 2001, 2002; Heim- 
pel et al. 2005; Ballot et al 2007; Kapyla et al. 2011; Gastine 
and Wicht 2012; Guerrero et al 2013; Hotta et al 2015). The 
gas planets Jupiter and Saturn (Porco et al. 2003, 2005), and the 
Sun (Thompson et al 2003) belong to this category. Significant 
prograde zonal flows might also be present in the Earth’s outer core 
near the equator (Gillet et al 2015). These zonal winds are excited 
due to the tilted columns of convection (seen in the vorticity plot) 
which arise due to the combined effect of Coriolis forces and the 
curvature of the outer boundaries in a spherical shell (Busse 1983; 
Zhang 1992; Busse 1994); however, also refer to Takehiro (2008) 
for a different interpretation. Such a convection pattern introduces 
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Non-magnetic cases with free-slip boundaries 



Non-magnetic cases with no-slip boundaries 



Dynamo cases with free-slip boundaries 
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Figure 2. Snapshots of axial vorticity in the equatorial plane and azimuthally averaged zonal flow in a meridional plane. The zonal flow velocity is expressed 
in terms of a Reynolds number. The shell is rotating anti-clockwise in the vorticity panels. The red and blue vorticity shades represent positive and negative 
values, respectively. In all the plots the color scale saturates below the extreme values to highlight fainter structures. Assuming a critical Ra of 7 x 10^ 
(Christensen and Aubert 2006) the supercriticality of convection is about 30 and 600 in the left and right column, respectively. 
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the so-called Reynolds stresses where a statistical correlation be¬ 
tween the azimuthal and the cylindrically-radial component of the 
flow gives rise to a transport of angular momentum towards the 
equatorial regions close to the outer boundaries. In passing we note 
that for rapidly rotating convection with high supercriticality (ap¬ 
proaching a ‘quasigeostrophic’ state) the Rhines scale determines 
the number of zones and amplitude of the zonal flow (Rhines 1975; 
Gastine et al 2014b). Here too, the curvature of the spherical shell 
enters via the so called “topographic /3 effect”. 

There are other mechanisms through which strong zonal 
flows can be excited. Systems with signiflcant latitudinal varia¬ 
tion of temperature/heat-flux can generate strong zonal flows via 
the ‘thermal-wind’ mechanism (e.g. Sreenivasan and Jones 2006; 
Miesch et al 2006). The so-called ‘compressional torque’ effect 
might also be an important ingredient in generating strong zonal 
flows in compressible convection (Glatzmaier and Gilman 1981, 
1982; Evonuk 2008; Glatzmaier et al 2009; Verhoeven and Stell- 
mach 2014); although, recent 3D anelastic simulations in spherical 
shells question the importance of this effect (Gastine et al 2014b). 

The vorticity structures are tilted in the opposite direction in 
panels Fig. 2b, where the Rayleigh number is strongly enhanced 
compared to the simulation shown in Fig. 2a. Moreover, the coher¬ 
ent structures with alternating cyclonic/anti-cyclonic vorticity seen 
in Fig. 2a are no-longer present. A rather large region of predomi¬ 
nantly anticyclonic vorticity (in blue) is visible in Fig. 2b whereas 
more fractured and chaotic structures All other regions of the shell. 
The zonal flow in Fig. 2b has changed sign and is now westward 
(retrograde) in low-latitude regions near the outer boundary. In this 
regime of vigorous convection the angular momentum is homoge¬ 
nized throughout the spherical shell (Aurnou et al 2007 ; Gastine 
et al 2013). The ice giant planets Uranus and Neptune might belo- 
ing to this category as they both show retrograde zonal flows (e.g. 
Fry et al 2012; Martin et al 2012). At this stage the Coriolis forces 
have a weaker influence on convection and the inertial effects be¬ 
come dominant. For even higher Ra, the convection eventually be¬ 
comes isotropic without any appreciable axisymmetric flow. This 
was indeed the case of the largest Ra cases with . 


3.1.2 Strong shear 


As discussed above zonal flows are a rather prominent feature in 
rotating convection in spherical shells. The competing influence of 
the buoyancy forces and the Coriolis forces dictate what would 
be the nature of the zonal-winds. Gilman (1977) introduced the 
so called convective Rossby number Roc — y/RoE^JPr which 
roughly quantifles the ratio of these two competing forces. It is in¬ 
structive to analyze the energy contained in the axisymmetric zonal 
flows in the setup discussed above. As a measure of the relative 
contribution of the zonal flows we use 





rptot 

-^kin 


(14) 


where and E]^^^ is the kinetic energy contained in the 

axisymmetric zonal flow component, in the non-zonal flow compo¬ 
nent, and in the total flow, respectively (Christensen 2002; Gastine 
and Wicht 2012). 

Variation of lu as a function of Roc is portrayed in Fig. 3 
for the non-magnetic cases with free-slip boundaries. As Roc in¬ 
creases for each Ekman number, an increase in m shows that frac¬ 
tion of energy in the zonal flow increases after the onset of the 
convection. After reaching a peak value m starts declining. This 



Figure 3. Variation of tu as a function of the convective Rossby number 
Roc for hydrodynamic cases with free-slip boundaries. The shaded zone 
highlights the zonal flow bistability region. 


may be attributed to the gradual disruption of the coherent con¬ 
vective columns when Roc increases, leading to reduced Reynolds 
stresses which transfer energy from small-scale convective motions 
to large-scale azimuthal flow (Christensen 2002; Gastine and Wicht 
2012). Around Roc ~ 1 a discontinuous change occurs in vj for 
E=l^~^ and 10“^, after which it starts to decline again. These 
trends in zonal flow strength have been observed in earlier simula¬ 
tions as well (Christensen 2002; Aurnou et al 2007; Gastine and 
Wicht 2012; Gastine et al 2013, 2014a). The sudden change in vj 
observed around Roc ~ 1 is due to the reversal of the prograde jet 
at lower Roc to a retrograde one at higher Roc. The zonal flow en¬ 
ergy in the retrograde-jet case exceeds that of the maximum energy 
in the prograde-jet case for E=10~^, consistent with Gastine et al 
(2013). However, the cases with E=10~^ shows the opposite trend. 
It is likely due to the large Ekman number in this setup, implying a 
more dominant role of viscosity. 

Convection in two dimensions with free-slip boundaries and 
one periodic direction is known for generating highly energetic 
‘winds’ or shear-flows in the direction of periodicity. Such winds 
are highly efficient at reducing the transport of heat perpendicular 
to its own direction (see e.g. Terry 2000; Garcia et al 2003; Go- 
luskin et al 2014). These systems are very much similar to our 
setup: the tendency of two dimensionality is introduced by the in¬ 
fluence of the Coriolis forces on the convection, the boundaries are 
similarly free-slip, and the azimuthal direction in spherical shell 
is akin to the periodic boundaries in Cartesian boxes. Therefore, 
it can be expected that shear might have a strong influence on the 
heat transfer in our simulations as well. The axisymmetric zonal- 
jets do not participate in radial heat transfer, and the more energy 
they contain the better they become at disrupting heat-transporting 
motions which might try to cross them. 


3.1.3 Zonal flow bistability 

With the qualitative idea of stronger zonal flows leading to less effi¬ 
cient heat transfer in mind, we turn back to a feature visible in Fig. 1 
which has not been discussed so far. For Ekman number of 10“^ 
and 10“^, regions where two different solutions coexist for cases 
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with free-slip boundaries, are highlighted with gray-shaded boxes 
in Fig. 1. In a recent study dedicated to stellar differential rotation 
Gastine et al (2014a) demonstrated the possibility of bistability 
in the nature of the differential rotation in rotating spherical shells 
with aspect ratio of 0.6. This result obtained in Boussinesq convec¬ 
tion was later confirmed by fully compressible simulations (Kapyla 
et al 2014); however, see Karak et al (2015). We investigate this 
feature for our simulation setup as well. Prograde zonal jet persists 
over a wide range of low and intermediate values of the Rayleigh 
number, but at a sufficiently high-enough Ra, the equatorial zonal- 
jet reverses direction and becomes retrograde. If we now take this 
state with retrograde zonal-jet and use it as an initial condition for a 
simulation with a smaller Ra the retrograde zonal-jet still persists 
and appears to be a stable feature. The reversal of the zonal-jet and 
the bistability phenomenon occurs when Roc ~ 1, as highlighted 
in Fig. 3 with a shaded region. In analogy with typical non-linear 
systems with hysteresis, there are two possible attractors for a given 
set of parameters in the bistable regime, and depending on the ini¬ 
tial conditions a simulation can either be captured by the prograde 
attractor or the retrograde one. This phenomenon only occurs in a 
rather limited Ra range. 

For the retrograde-jet solution has somewhat smaller 

energy in the zonal fiow, while, for E=10~^, the prograde-jet solu¬ 
tions carries less energy in the zonal flow. Therefore, it appears that 
more energy in the retrograde-jet solution for E=l{)~^ leads to a 
lower Nusselt number as compared to the corresponding prograde- 
jet solutions. This scenario is reversed for E=l^~^ simulations. 
Therefore, the difference in zonal flow energy in Fig. 3 qualitatively 
explains the different Nu 'm the gray-shaded boxes in Fig. 1. 

3.1.4 Latitude dependent Nusselt number 

The distribution of heat flow with latitude on the outer boundary 
provides interesting information on the flow regimes inside and 
outside the tangent cylinder (TC). In Fig. 4a we plot a latitude- 
dependent Nusselt number 

= (atr = r„) (15) 

on the outer boundary for some selected HD cases with E=l^~^ 
and free-slip boundaries (solid curves). At low Ra a peak in Nu(0) 
in the equatorial region (small latitudes) shows that the regions 
inside the TC remain rather quiescent for small convective driv¬ 
ing (Jones et al 2000; Dormy et al 2004). As Ra increases, the 
heat-flux at high latitudes increases rapidly once convection in¬ 
side the TC sets in. Remarkably, the heat-flux in the equatorial re¬ 
gions nearly saturates (compare the blue and the green solid curves 
in Fig. 4a). This is due to the severe quenching of radially out¬ 
ward flows by the strong zonal jets. In the equatorial regions, the 
zonal winds are nearly perpendicular to the direction of gravity, 
and, hence, the maximum impact is felt in this region. As we go 
to higher latitudes this effect decreases. Hence, we can conclude 
that heat-transfer efficiency near the equator is highly affected by 
zonal-jets. 

Figure 4a highlights another interesting aspect. The rise in 
heat transfer with increasing Rayleigh number proceeds much 
faster in the polar regions (once TC convection has started) than it 
does in the equatorial regions. Even though convection in the equa¬ 
torial regions starts at lower Ra, the polar Nusselt number quickly 
rises above the equatorial one at higher values of Ra. It qualita¬ 
tively implies that in a relation of type Nu oc Ra^ (dependence 
on other parameters is suppressed for the sake of simplicity) the 


exponent a will be higher in the polar regions as compared to the 
equatorial one. It may provide a plausible explanation why a re¬ 
ported in rotating box simulations, which are suited for the polar 
regions in a spherical shell, is generally higher (King et al 2012; 
Ecke and Niemela 2014; Cheng et al 2015) than the one reported 
in a similar study in spherical shells (King et al 2010). 


3.2 Models with no-slip boundaries 

Elow boundary conditions have a strong effect on the zonal flow 
strength and nature (Gilman 1978; Christensen et al 1999; Chris¬ 
tensen 2002; Aubert 2005; Busse and Simitev 2006; Schrinner et al 
2012; Dharmaraj et al 2014; Garcia et al 2014). Here we inves¬ 
tigate the effect of changing the mechanical boundary conditions 
from free-slip to no-slip on the nature of heat-transport. In Eig. 1 
we also plot Nu vs. Ra for non-magnetic simulations with no¬ 
slip boundaries (empty data points). Qualitatively both no-slip and 
free-slip data sets follow similar trends, in that there are 3 different 
regimes (onset, rotation dominated, and inertia dominated). There 
are two important differences which can be easily identified: in the 
onset and in the rotation dominated regime, the Nusselt number is 
higher in the simulations with no-slip boundaries, while it gradually 
becomes lower than the free-slip cases as Ra is further increased. 
Differences in Nu between cases with free-slip or no-slip bound¬ 
aries has been reported in earlier weakly-rotating RBC studies. In 
these setups a higher value of in the free-slip cases may be 
attributed to the absence of viscous boundary layers (Julien et al 
1996). 


3.2.1 Flow structure 

The axial vorticity and the zonal flow plots for the no-slip HD cases 
are shown in Eig. 2(c,d). A substantial change in both panels, as 
compared to their free-slip counterparts in Eig. 2(a,b), is immedi¬ 
ately noticeable. The much weaker zonal flows show more variation 
along the rotation axis in the no-slip cases. Due to the reduced zonal 
shear, the spiralling behavior of the vorticity structures is greatly 
reduced. Another point to notice is that the dichotomy of prograde 
and retrograde zonal jets on the outer boundary and near the equa¬ 
tor is much less apparent here. The reason for the quenched zonal 
flows can be attributed to the presence of Ekman boundary layers 
with no-slip boundaries (e.g. Jones 2007). In simulations with free- 
slip boundaries the saturation of the zonal flows is controlled by 
the weak bulk-friction due to viscosity, while, in the case of no-slip 
boundaries, the much larger friction in the Ekman boundary layers 
((X see Gillet and Jones (2006); Jones (2007)) inhibits vigor¬ 
ous zonal flows more efficiently (Christensen 2002; Aubert 2005; 
Jones 2007). 


3.2.2 Latitude dependent Nusselt number 

Eigure 4a also displays Nu(0) for no-slip HD simulations with 
dashed curves. Inside the TC no-slip HD cases resemble their free- 
slip counterparts. However, in the low-latitude regions outside the 
TC, the difference is rather substantial. As discussed before, the 
zonal flow is quenched by the no-slip boundaries, and, therefore, it 
does not impede the heat transfer as efficiently as in the free-slip 
case. As a result, the low-latitude regions in no-slip simulations re¬ 
main efficient at transporting heat. 
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Figure 4. Variation of the latitude-dependent Nusselt number Nu{6) (see Eq. 15) at the outer boundary for simulations with Ekman number of 10“^. Different 
colors in each panel represent different Rayleigh numbers: red for Ra=S x 10^, blue for Ra=2 x 10®, and green for i?a=6 x 10®. The trail of blue triangular 
points represent the bistable multipolar dynamo solution for the same control parameter as those used for the blue dashed curve (dipolar solution). The vertical 
broken gray-lines represent the latitudes where the tangent cylinder touches the outer boundary. Assuming a critical Ra of 10^ (Christensen and Aubert 2006) 
the supercriticality of convection is about 8, 20, 60 for the three Ra mentioned above. 


3.3 Relation between shear flow and heat transport 


To show the effect of zonal flow on the heat-transfer we compare 
in Fig. 5 two relevant quantities as a function of Roc. One is the 
ratio of the Nusselt number in the no-slip case, Nun, to the cor¬ 
responding value in the free-slip case, iVu/, at otherwise the same 
parameter values. The other, distinguished by black-rimmed sym¬ 
bols, is the parameter vjf for free-slip cases (our measure for the 
relative contribution of zonal flow to the kinetic energy). Note that 
zu for the no-slip cases remains close to unity due to quenched 
zonal flows. As can be clearly seen in Fig. 5 there is a tight cor¬ 
relation between the two quantities. Qualitatively, it implies that 
the larger the amount of energy contained in the zonal jets in the 
free-slip cases, the larger will be the impact on heat-transfer if we 
quench these jets by replacing the free-slip boundary condition by 
no-slip. A sudden increase in the iVu-ratio near Roc ~ 1.5 is in¬ 
duced by the increase in the energy content of the zonal flow due to 
the transition to the retrograde-jet regime in the free-slip cases. In 
the large Roc regime zu is expected to approach unity since the sys¬ 
tem will resemble an RBC type state where any axisymmetric flow 
will carry negligible energy as compared to the total kinetic energy. 
This is starting to happen for the E=10~^ cases. The A^u-ratio, 
however, will not be unity as long as the influence of the boundary 
layers is felt by convection; it should be smaller than unity for high 
Roc since free-slip cases transport more heat (Julien et al 1996). 

It should be noted that we have not analyzed the contribu¬ 
tion of the Ekman pumping mechanism, which is present in no¬ 
slip cases, on the heat-transfer efficiency gain (e.g. Stellmach et al 
2014). However, noting that Nu was very much similar in both 
free-slip and no-slip HD simulations inside the TC we may specu¬ 
late that this effect is smaller at our control parameters. 

To summarize, in the rotation dominated regime (Roc < 1), 
the simulations with no-slip boundaries transfer more heat because 
of their much weaker zonal flows. In the inertia/buoyancy domi¬ 
nated regime (Roc > 1), however, the cases with free-slip boud- 
naries transfer more heat due to the lack of viscous boundary layers. 
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Figure 5. On the left axis, ratio of the Nusselt number Nu for no-slip (sub¬ 
script ‘n’) and free-slip (subscript ‘f’) cases, shown by data points without 
rim, and, on the right axis, zuf (black rimmed data points) for free-slip 
cases as a function of the convective Rossby number Roc. The data set is 
for non-magnetic simulations with aspect ratio of 0.35. 


3.4 Effect of shell thickness 

The dominance of zonal flows appears to be somewhat enhanced 
in spherical shells with larger aspect ratios and free-slip bound¬ 
aries (e.g. Aurnou et al 2008; Gastine et al 2013). This would im¬ 
ply that the difference in the Nusselt number between free-slip and 
no-slip simulations would be larger in spherical shells with thinner 
convection zones. To test this idea we have run simulations with an 
aspect ratio of 0.6. In Fig. 6 we plot the Nu ratio and free-slip 'uj 
for this new data set. As expected the enhancement of convective 
efficiency in the no-slip case compared to the free-slip case is even 
higher in this setup, increasing about three fold in some cases. 
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Figure 6. On the left axis, ratio of the Nusselt number Nu for no-slip (sub¬ 
script ‘n’) and free-slip (subscript ‘f’) cases, and, on the right axis, vuf 
(black rimmed data points) for free-slip cases as a function of the convec¬ 
tive Rossby number Roc. The data set is for non-magnetic simulations with 
aspect ratio of 0.6. Note the additional data points at E=3 x in gray. 




Figure 7. Spectrum of kinetic energy and magnetic energy as a function of 
the spherical harmonic degree i. The data is for a snapshot of the simulation. 
Panel a and b correspond to Fig. 2e and 2f, respectively. 


4 DYNAMO SIMULATIONS 

In the previous section we showed that the presence of zonal flow is 
an important factor for the heat transfer efficiency in rotating con¬ 
vection and that the relative magnitude of the zonal jets depends on 
the mechanical boundary conditions. The magnetic held generated 
by a dynamo process can also, if it is strong enough, suppress zonal 
flows (see e.g. Aubert 2005; Miesch 2005; Browning 2008; Yadav 
et al 2013b; Heimpel and Gomez Perez 2011; Duarte et al 2013; 
Karak et al 2015). To explore this aspect we run self-consistent 
dynamo simulations for both kinds of mechanical boundary condi¬ 
tions. 

4.1 Dynamos with free-slip boundaries 

We begin by analyzing dynamo simulations with free-slip boundary 
conditions. Comparing these simulations with the corresponding 
HD simulations with free-slip boundaries will help us to pin down 
the change in heat-transfer efficiency due to the interaction of the 
magnetic held with the zonal-flows. 

4.1.1 Flow structure 

In Fig. 2(e,f) we display the axial vorticity and the zonal flow pro¬ 
files for dynamo simulations with free-slip boundaries. The zonal 
flow strength is reduced in Fig. 2e as compared to its HD counter¬ 
part in Fig. 2a. The reduction in zonal flow is less pronounced for 
the high-Ra case in Fig. 2f. As a result of the reduced zonal shear, 
the vorticity structures seen in Fig. 2e are more radially oriented in 
the equatorial plane. Both of these effects are qualitatively similar 
to what we already observed when we changed from free-slip to 
no-slip boundaries in HD simulations. 

A feature similar to the broad anticyclonic vorticity region 
seen in Fig. 2b is visible again in Fig. 2f. Soderlund et al (2013) 
have also reported such vorticity structures in longitudes (see their 
‘thick-dynamo’ case). The magnetic held is highly efficient at 
quenching the shear flow only when it is of sufficiently large scale 


and has energy comparable to or larger than the kinetic energy. In 
the dynamo simulations considered here the total magnetic energy 
(ME) is close to the total kinetic energy for cases with Roc < 1. 
This is due to the characteristic columnar convection which sup¬ 
ports the generation of a large-scale dynamo. In this regime, the 
magnetic morphology could be either of geodynamo-type where 
the axial dipolar component carries the largest share of the total 
ME (Kageyama and Sato 1997; Olson et al 1999; Aubert et al 
2008), or it could be multipolar where higher-order held com¬ 
ponents dominate over the axial dipolar one (Busse and Simitev 
2006; Christensen and Aubert 2006; Brown et al 2010; Kapyla 
et al 2012; Gastine et al 2012). Typically, in both kinds of dy¬ 
namo solutions, the magnetic held at large scales remains compara¬ 
ble to the kinetic energy. Eigure 7a shows the comparable strength 
of the kinetic and the magnetic energy at different spherical har¬ 
monic degrees i for the multipolar dynamo case shown in Eig. 2e. 
Systematic studies of the geodynamo have revealed that the mag¬ 
netic energy becomes increasingly dominant as the Ekman number 
is decreased and can be an order of magnitude larger than the ki¬ 
netic energy (Christensen and Aubert 2006; Takahashi and Shimizu 

2012) . Therefore, for lower E, the ratio ME/KE is likely to increase 
even further, enhancing the influence of the magnetic held on the 
flow. On the other hand, when Roc is greater than unity, the he¬ 
lical columns are no longer dominant and plume-like convection 
generates a small-scale dynamo powered by the local chaotic mo¬ 
tions. The spectrum shown in Eig. 7b, which corresponds to the 
case given in Fig. 2f, highlights the much reduced magnetic energy 
at low i or large length scales. Therefore, we can expect that at 
small Roc, the back-reaction of the magnetic held on the flow will 
be rather strong, but as Roc becomes larger, the dynamo becomes 
less efficient and the influence of the magnetic held would become 
much weaker. This explains why the prograde zonal-jets are highly 
quenched for low Roc dynamos while, for larger Roc, the retro¬ 
grade zonal jets are not dramatically influenced (Soderlund et al 

2013) . 
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4.1.2 Latitude dependent Nusselt number 

The variation of heat flux with latitude for free-slip dynamo cases 
is shown in Fig. 4b with dashed curves, along with the correspond¬ 
ing profiles for free-slip hydrodynamic cases. We again see that 
due to the quenching of the zonal flow by the magnetic held the 
heat transport at low latitudes is enhanced. However, unlike the HD 
simulations with no-slip boundaries, a significant enhancement of 
Nu even inside the TC can be observed in Fig. 4b. The Nusselt 
number enhancement inside the TC only works for high enough 
Ra since it is not present for lowest Ra in Fig. 4b. 

The dynamo simulation for which Nu\^ enhanced inside the 
TC produces a dipole dominant magnetic held. This configuration 
has almost vertical poloidal magnetic held threading through the 
region inside the TC. We note that dynamo simulations with free- 
slip boundaries exhibit an interesting properly called “bistability” 
for a limited range of control parameters (Simitev and Busse 2009; 
Schrinner et al. 2012; Gastine et al. 2012). Such setups can sustain 
either a dipole dominant magnetic held or a multipolar magnetic 
held depending on the configuration and strength of the magnetic 
held used as initial condition. We exploit this property and gener¬ 
ate a multipolar dynamo solution for the same control parameters 
as those used in generating the dipolar solution at Ra — 2x 10®. In 
Fig. 4b we show the Nu(0) profile for this multipolar dynamo so¬ 
lution using a trail of filled triangles. In this case Nui^ reduced in¬ 
side the TC as oppose to an inhancement in the dipolar case shown 
with dashed curve. Close to the equator, the profiles are similar for 
dipolar and multipolar solutions. Therefore, for some simulations, 
the Nu enhancement is due to two reasons: zonal flow quenching 
at low-latitudes and significant modifications of convection at high 
latitudes by the magnetic held. 

4.2 Relation between shear flow and heat transport 

Similar to Fig. 5 we compare the Nu ratio for free-slip HD and 
free-slip dynamo simulations and vj for free-slip HD simulations 
in Fig. 8. Once again a good correlation exists between the two 
quantities, re-enforcing the idea that quenching of zonal jets is the 
primary cause of the heat-transfer efficiency enhancement in the 
dynamo cases compared to the non-magnetic cases. 

For a Rossby number larger than unity the Nusselt number 
of dynamo simulations is somewhat smaller than that of the cor¬ 
responding HD simulations. This indicates that the small-scale dy¬ 
namo operating in the high Roc regime is introducing the magneto¬ 
quenching effect where the presence of the magnetic held makes 
the convection less efficient at transferring heat. Since the magnetic 
held itself is weak in this regime, the associated decrease in iVu is 
also rather small. 

It is possible that the magnetic quenching effect is present even 
for Roc < 1. In this case, even though the magnetic held may 
quench convection to some extent, the increase in Nu due to the 
reduced zonal flow is much larger than the reduction 'm Nu due 
to magnetic quenching. For larger Roc, on the other hand, since 
the effect of zonal flow is much weaker, the reduction in by 
magnetic quenching prevails. 

4.3 Dynamos in deep convective layers 

We also compare the heat-transfer efficiency in HD and dynamo 
simulations in spherical shells with very deep convection zones 
(shells with aspect ratio of 0.2). For this purpose we analyse dy¬ 
namo cases from Gastine et al. (2012), supplemented by additional 



Ro^ 


Figure 8. On the left axis, ratio of the Nusselt number Nu for free-slip dy¬ 
namo (subscript ‘d’) and free-slip HD (subscript ‘h’) cases, and, on the right 
axis, vjf (black rimmed data points) for free-slip HD cases as a function of 
the convective Rossby number Roc. The A/'n-ratio for the bistable dipolar 
dynamos in the range 0.2 < Roc < 0.3 is also plotted using small empty 
triangles. The data set is for simulations with aspect ratio of 0.35. 


similar hydrodynamic simulations. These cases are not completely 
analogous to our basic setup: the inner boundary is no-slip while 
the outer boundary is free-slip. Since the non-convective inner core 
is rather small its boundary condition is unlikely to strongly affect 
the dynamics of the global convective layer. 

We plot Nu ratio and m in Fig. 9. The contribution of zonal 
flow is reduced in this thick-shell configuration as compared to the 
earlier thinner spherical shells. For instance, the parameter vjh, 
defining the relative stregnth of the zonal flow, reaches a maxi¬ 
mum value of about 25 at E=10~^ in the HD simulations with 
?7=0.2 while it reached about 50 in the HD cases with ?7=0.35 and 
E=10~^ (Fig. 8). It shows that the zonal flow strength is weaker 
in the thick shell cases. For E=10~^ cases, there is no correlation 
between the Nu ratio and Wh- For E=l^~^ a correlation exists 
but only for large enough Roc, and, at E=10~^, a tight correla¬ 
tion for low Roc and a rough correlation at higher Roc appears. 
Therefore, as we decrease the Ekman number, the energy content 
of the zonal flow increases and we start to see a better correlation 
between Nu ratio and vjh- Nonetheless, comparing this data-set 
with that for p = 0.35 (Fig. 8) indicates that spherical shells with 
deep convective layers promote weaker zonal flows as compared to 
thinner spherical shells. The weaker enhancement in A^u due to the 
presence of a dynamo-generated magnetic held follows as a conse¬ 
quence. 


4.4 Dynamos with no-slip boundaries 

In the earlier sections we saw that either the no-slip mechanical 
boundary conditions or the presence of a dynamo-generated mag¬ 
netic held enhance the convective efficiency by quenching the zonal 
shear. The magnetic held, however, is special in that it not only af¬ 
fects the zonal shear but also changes the structure of the convective 
flow responsible for transferring heat. In this sub-section we in¬ 
vestigate dynamo simulations which have no-slip flow boundaries. 
Comparing the no-slip HD simulations with these dynamo Simula- 
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Figure 9. On the left axis, ratio of the Nusselt number Nu for free-slip 
dynamo (subscript ‘d’) and free-slip HD (subscript ‘h’) cases, and, on the 
right axis, vjf (black rimmed data points) for free-slip HD cases as a func¬ 
tion of the convective Rossby number Roc. The data set is for simulations 
with aspect ratio of 0.2. 


tions will allow us to minimize the effects of different zonal flow 
amplitudes on the Nusselt number. 


4.4.1 Flow structure 

The axial vorticity plots and the zonal flow profiles are shown in 
Fig. 2(g,h) for these new dynamo cases. The zonal flows are highly 
quenched, and, correspondingly, the shearing in the axial vorticity 
structures is much weaker. The zonal flow also varies substantially 
along the rotation axis, usually observed in dynamos with dipole- 
dominant magnetic held configurations (Aubert 2005; Schrinner 
et al. 2012; Yadav et al. 2013b). In the high Ra regime (small- 
scale dynamo cases), the zonal flow profile is similar to the cor¬ 
responding hydrodynamic case in Fig. 2d. Once again it illustrates 
that small-scale dynamos have a much smaller influence on the flow 
as compared to the dipolar or the multipolar dynamos at smaller 
Roc (at least in our simulations). 


4.4.2 Latitude dependent Nusselt number 

The latitudinal heat-flow profiles for no-slip dynamo cases are 
shown by dashed curves in Fig. 4c, along with the correspond¬ 
ing profiles for no-slip HD cases. Near the equator both HD and 
dynamo cases behave somewhat similarly. This is due to the fact 
that the zonal jets, responsible for controlling the heat-transfer at 
low-latitudes in the free-slip cases, are highly quenched in both 
setups. Inside the TC, however, different behaviour exist at dif¬ 
ferent Rayleigh numbers: similar heat flow for Ra=^ x 10^ (red 
curves), increased heat flow in the dynamo case at Ra=2 x 10® 
(blue curves), and mixed behavior for Ra=Q x 10® (green curves). 
The sharp increase close to the poles in Fig. 4c appears to be due 
to the presence of large-scale helical polar vortices with enhanced 
radial and horizontal flows as compared to the corresponding HD 
case. Such vortices are preferably excited when a strong dipolar 
held is present (Sreenivasan and Jones 2005). 


4.4.3 Relation between field strength and heat transport 


In both HD as well as dynamo simulations with no-slip boundaries 
the energy in the zonal flow is rather small. Therefore, the zonal 
flows do not play an important role in determining the heat-transfer 
efficiency in these setups. In Fig. 10 we compare the ratio of Nu 
for dynamo and HD cases as a function of the mean magnetic held 
strength {B) which is quantified by the classical non-dimensional 
parameter called the Elsasser number 


PfiXQ, 


(16) 


This figure shows two distinct features. For the multipolar dynamos 
(empty symbols), the magnetic held only marginally quenches the 
heat transport. For the dipolar dynamos at E=10~^, however, a 
peak is present at A 3 where the efficiency of convection in 
transporting heat is enhanced by 30%. At a higher Ekman number 
of 10~^ the enhancement of the Nusselt number is smaller and re¬ 
verses for A > 3. Note that the dashed line for E=10~^ data points 
follows the increasing Roc trend. When the dynamo changes be¬ 
havior from dipolar to multipolar the held strength decreases even- 
though Roc is increased. Therefore, the blue dashed-line does not 
follow a strictly increasing trend. 

The shape of the green curve (for E=10~^) is quite similar 
to that reported by King and Aurnou (2015). In this study liquid 
Gallium (Pr ~ 0.025 and Pm ~ 10“®) was used as the work¬ 
ing fluid and rotating convection (2 x 10“® ^ E ^ 2 x 10“^) 
was investigated in a cylindrical container with/without a verti¬ 
cally imposed magnetic field. This study reported maximal con¬ 
vective efficiency (about 50% higher Nu in magnetoconvection) at 
A ^ 2 and E ^ 10“®. The authors interpret the A^u-enhancement 
as a signature of ‘magnetostrophic’ regime of convection (Chan¬ 
drasekhar 1954; Stevenson 1979). Rotating convection in cylindri- 
cal/Cartesian setups with vertically imposed magnetic held can be 
compared with the polar-region convection in rotating spherical- 
shells generating a dipole-dominant magnetic held. 

The similarity of Nu enhancement in our E=10~^ cases with 
the experimental results by King and Aurnou (2015) suggests that 
magnetic held is playing a dominant role in governing the convec¬ 
tion behavior in the no-slip dynamo cases, especially at high lati¬ 
tudes (Eig. 4c). In this context, it is interesting to note a recent study 
by Teed et al. (2015) who analysed torsional oscillations in a mag¬ 
netoconvection setup in rotating spherical shells and And a possible 
transition to magnetostrophic regime also at E ^ 10“®. This hints 
that such Ekman number simulations might be entering the magne¬ 
tostrophic convection regime. However, without a detailed analysis 
of various forces in our simulations we can not confirm the exis¬ 
tence of such a regime. We defer this exercise to our forthcoming 
study where lower Ekman number simulations will be analysed. 


5 SUMMARY AND OUTLOOK 

In this study we sought to shed some light on the question of what 
factors affect the heat-transfer efficiency of convection in rotating 
spherical shells as it occurs in fluid regions of planets and stars. We 
compare four setups, hydrodynamic simulations with either free- 
slip or no-slip flow boundary conditions and the corresponding self- 
consistent dynamo cases. The key features can be summarised as 
follows: 

• Hydrodynamic simulations with free-slip (stress-free) bound¬ 
aries generate the strongest zonal jets. These axially-aligned zonal 
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Figure 10. Variation of the ratio of Nu for dynamo (subscript ‘d’) and 
hydrodynamic (subscript ‘h’) cases with no-slip boundaries as a function of 
the Elsasser number A. Filled/unfilled symbols represent dipolar/multipolar 
dynamos. Only dynamo simulations with Pm=l are shown. 


flows disrupt and impede the heat-transporting motions. Their ef¬ 
fect is most efficient in the equatorial regions where they severely 
quench the efficiency of convection, whereas, in the polar regions, 
the effect of the zonal flows on the heat-transport efficiency was 
substantially reduced. 

• Changing the free-slip boundary condition to no-slip strongly 
reduces the energy of the zonal flows. Consequently, the Nusselt 
number increases substantially, by as much as a factor of 3 in some 
of our simulations with a thin convecting shell. The heat transport is 
most strongly enhanced at low latitudes. However, in the inertia or 
buoyancy dominated regime at large Rossby numbers, the Nusselt 
number is slightly higher with free-slip boundaries. 

• The presence of a dynamo-generated magnetic held in the case 
of free-slip boundaries also efficiently quenches the zonal flows and 
enhances the heat flow at low latitudes. The heat-flux in the polar 
regions is affected by the magnetic held as well. Dipolar dynamos 
can substantially increase the polar heat-flux, whereas multipolar 
dynamos can quench it. Again, this holds in the rotation-dominated 
regime at low Rossby numbers. 

• In case of no-slip boundaries, when zonal flows play a mi¬ 
nor role also without a magnetic held, the presence of a dynamo¬ 
generated held can enhance the Nusselt number signiflcantly. Con¬ 
ditions for this to happen are that the Ekman number is low enough 
{E ^ 10“^) and the magnetic held is dipole-dominated and has a 
suitable strength, that is, the Elsasser number A is of order one. At 
E — 10“^ the optimum is reached for A 3. 

• Our results demonstrate the importance of self-excited zonal 
flows and magnetic flelds for the heat transfer efficiency of rotat¬ 
ing convection, as well as the different influence that they have in¬ 
side and outside the tangent cylinder. This suggests that studies of 
heat transfer in Cartesian layers (for instance. Barker et al 2014) 
or cylindrical geometry (for instance, Cheng et al 2015) may miss 
many crucial effects that are important in planetary or stellar appli¬ 
cations. 

There are many fronts where our study can be extended. We 
have not assessed the influence of the (magnetic) Prandtl number 
on our results. Increasing and decreasing the fluid Prandtl num¬ 
ber has a direct effect on the strength of the zonal flows in ro¬ 


tating shells, with zonal flows being stronger for Pr < 1 (e.g. 
Zhang 1992; Aubert et al 2001; Christensen 2002; Simitev and 
Busse 2005; Gillet et al 2007; Calkins et al 2012). The mag¬ 
netic Prandtl number might also be important since it affects the 
held strength and morphology of the dynamo-generated magnetic 
held (e.g. Christensen and Aubert 2006; Schrinner et al 2012; Ya¬ 
dav et al 2013b,a; Raynaud et al 2014). Compressibility, which is 
very important for convection in giant planets and stars, is another 
aspect which likely influences the convective efficiency by chang¬ 
ing the nature of the zonal flows (Gastine and Wicht 2012). Explo¬ 
ration of these aspects and testing the approach towards a magne- 
tostrophic flow regime by simulations at even lower values of the 
Ekman number are essential steps towards a better understanding 
of rotating and magnetic convection in natural objects. 
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